データセット

library(tidyverse)
path <- "~/R/project/R_tutoring/Git/R_tutoring/Data/kadai01.csv"

df_raw <- read_csv(path, locale = locale(encoding = "UTF-8"))

kadai01 <- df_raw %>% 
  mutate(pref = c("Hokkaido", "Aomori", "Iwate", "Miyagi", "Akita", "Yamagata", 
   "Fukushima", "Ibaraki", "Tochigi", "Gunma", "Saitama", "Chiba", 
   "Tokyo", "Kanagawa", "Niigata", "Toyama", "Ishikawa", "Fukui", 
   "Yamanashi", "Nagano", "Gifu", "Shizuoka", "Aichi", "Mie", "Shiga", 
   "Kyoto", "Osaka", "Hyogo", "Nara", "Wakayama", "Tottori", "Shimane", 
   "Okayama", "Hiroshima", "Yamaguchi", "Tokushima", "Kagawa", "Ehime", 
   "Kochi", "Fukuoka", "Saga", "Nagasaki", "Kumamoto", "Oita", "Miyazaki", 
   "Kagoshima", "Okinawa")
  )


# show
DT::datatable(kadai01, 
              rownames = FALSE, 
              extensions = 'Buttons',
              options = list(autoWidth = FALSE,
                             pageLength = 5,
                             dom = 'Bfrtip',
                             buttons = list("csv"),
                             scrollX = TRUE,
                             scrollCollapse = TRUE),
              class = 'cell-border stripe',
              caption = "source: ")

データセットの呼び出しを行う

データの読み込み(インポート)

  1. データの読み込み(※ディレクトリについて) ここでは,主にCSV(Comma Separated Values)ファイルの読み込みについて説明する。 (注)ディレクトリについて ・Windowsではメニュー「ファイル」→「ディレクトリの変更」で指定できる ・Macではメニュー「その他」→「作業ディレクトリの変更」 実際に,それぞれの環境の中でデータを読み込むときは,ディレクトリの理解が不可欠になる。 例えば,JドライブのRフォルダにある”sheet1.csv”を読み込むときは次のようにする (学院のPC教室ではJドライブに資料を格納しておく)。
# データを読み込む
kadai01 <- read.csv("J:/R/kadai01.csv")

1. データの提示

# show the first 6 lines of `kadai01`
knitr::kable(head(kadai01))
pref pop epc
Hokkaido 5384 11070
Aomori 1309 2720
Iwate 1280 2769
Miyagi 2334 4819
Akita 1023 2177
Yamagata 1123 2385

2.基本統計量の計算

平均値・最大値・最小値・中央値

summary(kadai01)
##      pref                pop             epc       
##  Length:47          Min.   :  574   Min.   : 1360  
##  Class :character   1st Qu.: 1114   1st Qu.: 2543  
##  Mode  :character   Median : 1649   Median : 3556  
##                     Mean   : 2704   Mean   : 5678  
##                     3rd Qu.: 2728   3rd Qu.: 5838  
##                     Max.   :13514   Max.   :28097

不偏分散

lapply(kadai01[,2:3], var)
## $pop
## [1] 7444685
## 
## $epc
## [1] 28866161

不偏標準偏差

lapply(kadai01[,2:3], sd)
## $pop
## [1] 2728.495
## 
## $epc
## [1] 5372.724

3.相関係数の計算

cor(kadai01[,2], kadai01[,3])
##           epc
## pop 0.9974586

4. 線形モデルの生成

kadai.lm <- lm(epc ~ pop, data = kadai01)

回帰直線を方程式を求める 

coef(kadai.lm)
## (Intercept)         pop 
##  365.697315    1.964112

上記の結果から、今回の回帰係数は、切片が365.6973149で、傾きが1.9641119であることがわかる。これを線形の式に直すと、以下の通りになる。

\[ epc = 365.6973149 + 1.9641119 \times pop \]

5. 散布図および回帰直線の描画

plot(kadai01[,2:3],
     xlab = "Tthe Number of Population (thousands)",
     ylab = "Electricity Consumption (million kWh")
abline(kadai.lm)

6.残差(テキストp.83)の絶対値が比較的大きい都道府県に注目し、自身の考察を書きなさい。

残差を求める

kadai01_result <- data.frame(kadai01, 
                      resid = resid(kadai.lm),
                      resid_abs = abs(resid(kadai.lm))
                      )

並び替える

head(kadai01_result[sort(kadai01_result$resid_abs, decreasing = F, index = T)[[2]],], n = 10)

# or
kadai01_result %>% 
  arrange(desc(resid_abs)) %>% 
  head(n = 10) %>% 
  knitr::kable()
kadai01_result %>% 
  arrange(desc(resid)) %>% 
  head(n = 10) %>% 
  knitr::kable()
pref pop epc resid resid_abs
13 Tokyo 13514 28097 1188.2946 1188.2946
34 Hiroshima 2845 6592 638.4044 638.4044
17 Ishikawa 1154 3159 526.7176 526.7176
16 Toyama 1067 2918 456.5953 456.5953
40 Fukuoka 5103 10825 436.4397 436.4397
33 Okayama 1922 4572 431.2796 431.2796
28 Hyogo 5537 11626 385.0152 385.0152
18 Fukui 787 2217 305.5466 305.5466
20 Nagano 2100 4710 219.6677 219.6677
26 Kyoto 2610 5703 210.9707 210.9707

7.自己ワーク

環境省が公表する「環境統計集(平成29年度版)」 のデータを用いて、グラフを1つ作成し、考察しなさい。枠内には、①使用データとそのURL、②プログラムソース、③グラフ、④考察の4点を記載しなさい。なお、データ整形は、必ずしもRを使用しなくても良い。必要があれば、データ整形に関する説明も、①使用データの部分に記載しなさい。

ファイル提出

コード等をwordに貼って、ファイルを提出する

APPENDIX

summry of kadai.lm

summary(kadai.lm)
## 
## Call:
## lm(formula = epc ~ pop, data = kadai01)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1318.15  -129.32   -35.57   174.31  1188.29 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 365.69731   79.91473   4.576 3.72e-05 ***
## pop           1.96411    0.02091  93.913  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 387 on 45 degrees of freedom
## Multiple R-squared:  0.9949, Adjusted R-squared:  0.9948 
## F-statistic:  8820 on 1 and 45 DF,  p-value: < 2.2e-16

予測値と残差

残渣(residuals)は、実際の値と予測値の差分で計算される。

以下の図だと、予測値の青の線と実際の値である赤の線の間が残差である。

# load packages
library(modelr)
library(scales)
library(plotly)

# define the data for visualization
kadai01_plot <- kadai01 %>%
  add_predictions(kadai.lm) %>% 
  add_residuals(kadai.lm) %>% 
  `names<-`(value = c("pref", "pop", "epc", "predicted_epc", "residuals")) %>% 
  pivot_longer(cols = c(-pref, -pop), names_to = "vars", values_to = "value")



# raw v.s. expected value
g1 <- 
  kadai01_plot %>% 
  filter(vars != "residuals") %>% 
  ggplot(aes(pop, value, color = vars, label = pref)) +
  geom_point() +
  geom_line() +
  scale_x_continuous(labels = comma) +
  scale_y_continuous(labels = comma) +
  labs(title = "The scatter plot of the predicted value and raw data based on `kadai01`",
       x = "# of population",
       y = "Electricity",
       color = "Name") +
  theme_minimal()

ggplotly(g1)

実際の値である人口の変化によって生じる残差の値を、グラフ化したのが以下の図になる。

# residual  

library(ggrepel)

kadai01_plot %>% 
  filter(vars == "residuals") %>% 
  ggplot(aes(pop, value, label = pref)) +
  geom_hline(yintercept =  0, colour = "black", linetype = "dashed") +
  geom_point() +
  geom_line() +
  geom_label_repel(nudge_x = TRUE, nudge_y = TRUE, check_overlap = TRUE) +
    labs(title = "The scatter plot of the residual",
       x = "Pop",
       y = "Residuals") +
  theme_minimal()

日本地図

pacman::p_load("NipponMap", "sf", "ggthemes")

map <- read_sf(system.file("shapes/jpn.shp", package = "NipponMap")[1],
               crs = "+proj=longlat +datum=WGS84")

map2 <- 
  kadai01_plot %>%
  pivot_wider(id_cols = c(pref, pop), names_from = vars, values_from = value) %>% 
  rename(name = pref) %>%
  left_join(map, by = "name") %>% 
  st_sf()
  


g3 <- 
  ggplot(map2, aes(fill = residuals, label = name)) + 
  geom_sf(size = .1) +
  labs(title = "The Residuals Accross Japanese Prefectures") +
  scale_fill_gradient2(
    name = "Residuals",
    labels = scales::number_format(big.mark = ","),
    low = "#15C6E3",
    mid = "white",
    high = "#EB7411",
    midpoint = 0) +
  theme_pander()

g3

# ggplotly(g3) %>%
#   widgetframe::frameWidget()